#include "plant_1DSOS_RDM_input_1Dsin.h"


// --------------------------------------------------
plant_1DSOS_RDM_input_1Dsin::plant_1DSOS_RDM_input_1Dsin(double idt, int iLx,int iMode,double iE_D,double iE_S,double iE_N,bool iCAL_MODE,int iAgg_Length):
  Lx(iLx),
  Lx_agg(iLx/iAgg_Length),
  E_D(iE_D),
  E_S(iE_S),
  E_N(iE_N),
  mode(iMode),
  CAL_MODE(iCAL_MODE),
  Agg_Length(iAgg_Length)
{  
  // cout << "plant_1DSOS_RDM_input_1Dsin::plant_1DSOS_RDM_input_1Dsin()\n";
  
  // Options -----------------------------------
  dt = idt;
  N_CL	= NUMBER_NEIGHBOR;
  assert(iMode<=Lx_agg/2+1);

  // BASIC BLOCK ------------------------------
  // Input#0 = W (layer/s)
  // Input#1 = A (layer/s)
  // Input#2 = Freq
  // Input#3 = T (K)
  int idim[4] = {1,1,1,1};
  // Output#0 = h (layer/s)
  // Output#1 = r2 (layer^2)
  // Output#2 = m2
  // Output#3 = a (site)
  // Output#4~#5 = values of eigenfunction phi_1 ~ phi_2
  int odim[6];
  odim[0] = 1;
  odim[1] = 1;
  odim[2] = 1;
  odim[3] = 1;
  odim[4] = (mode+1); // The mode with sin
  odim[5] = (mode+1); // The mode with cos
  block::initialize(4,idim,6,odim,idt);

  // States ------------------------------------
  // Initialize random number generator --------
  gsl_rng_env_setup();
  gsl_T = gsl_rng_default;
  gsl_r = gsl_rng_alloc(gsl_T);
  gsl_rng_set(gsl_r,time_seed());

  work = gsl_fft_real_workspace_alloc(Lx_agg);
  real = gsl_fft_real_wavetable_alloc(Lx_agg);

  migration_events  = new MCArray(Lx,gsl_r);
  desorption_events = new MCArray(Lx,gsl_r);
  adsorption_events = new MCArray(Lx,gsl_r);
  ads_mig_des       = new MCArray(3,gsl_r);
  height            = new int[Lx];
  SQRT_ONE_OVER_LX = sqrt(1/double(Lx));
  SQRT_TWO_OVER_LX = sqrt(2/double(Lx));
  SQRT_ONE_OVER_LX_AGG = sqrt(1.0/double(Lx_agg));
  SQRT_TWO_OVER_LX_AGG = sqrt(2.0/double(Lx_agg));

  // Variables fro aggregationi
  height_agg        = new double[Lx_agg];

  reset();

  // Save a log --------------------------------
  ofstream flog("plant_1DSOS_RDM_input_1Dsin.log");
  flog << "Lx         = " << Lx  << "\n";
  flog << "E_D        = " << E_D << "\n";
  flog << "E_S        = " << E_S << "\n";
  flog << "E_N        = " << E_N << "\n";
  flog << "Agg_Length = " << Agg_Length << "\n";
  flog << "Lx_agg     = " << Lx_agg << "\n";
  flog.close();
  // cout << "END OF plant_1DSOS_RDM_input_1Dsin::plant_1DSOS_RDM_input_1Dsin()\n";
}
// --------------------------------------------------
plant_1DSOS_RDM_input_1Dsin::~plant_1DSOS_RDM_input_1Dsin()
{
  // cout << "Begin of ~plant_1DSOS_RDM_input_1Dsin()\n";
  delete migration_events;
  delete desorption_events;
  delete adsorption_events;
  delete ads_mig_des;
  delete []height;

  // Aggregation variables
  delete []height_agg;

  gsl_rng_free(gsl_r);
  gsl_fft_real_wavetable_free(real);
  gsl_fft_real_workspace_free(work);
  // cout << "End of ~plant_1DSOS_RDM_input_1Dsin()\n";
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::reset(){
  block::time = 0.0;

  // States
  for(int x=0;x<Lx;x++){
    height[x] = 0;
  }
  for(int x_agg=0;x_agg<Lx_agg;x_agg++){
    height_agg[x_agg] = 0.0;
  }
  
  // Inputs
  T	= 480;
  set_temperature(T);
  W	= 1.0;
  Freq	= 5;
  A	= 0.0;
  set_ads_rate(W,Freq,A);
  // Outputs
  h  = 0.0;
  r2 = 0.0;
  m2 = 0.0;
  a  = Lx;
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::update(double sysTime){
  //ofstream f_phi_definition("phi_definition.dat");
  //ofstream f_phi_FFT("phi_FFT.dat");

  if(time>sysTime)
    return;

  // Update states
  double iW	  = inputs[0][0];
  double iA	  = inputs[1][0];
  double iFreq    = inputs[2][0];
  double iT	  = inputs[3][0];

  if(iT!=T){
    set_temperature(iT);
  }
    
  if(iW!=W || iFreq!=Freq || iA != A){
      set_ads_rate(iW,iFreq,iA);
  }

  while(time<=sysTime){
    /*
    rate_ttl = ads_ttl+migr_ttl+desorp_ttl;
    double rdm = RandRand()*rate_ttl;
    if(rdm<ads_ttl)
      adsorption();
    else if(rdm<=migr_ttl+ads_ttl)
      migration();
    else
      desorption();
    */
    
    int event_id = ads_mig_des->PickAnEvent();
    switch(event_id){
      case 0:
        adsorption();break;
      case 1:
        migration();break;
      case 2:
        desorption();break;
    }
    time += OneMCTime(ads_mig_des->GetTotal());
  }

  // Calculate aggregation
  // Update the aggregation height
  for(int x_agg=0;x_agg<Lx_agg;x_agg++){
    height_agg[x_agg] = 0.0;
    for(int x_section=0;x_section<Agg_Length;x_section++){
      height_agg[x_agg] += height[x_agg*Agg_Length+x_section];
    }
    height_agg[x_agg] = height_agg[x_agg]/Agg_Length;
  }

  outputs[0][0] = getH();
  outputs[1][0] = getR2();
  outputs[2][0] = getM2();
  outputs[3][0] = 0.0; // FIXME, always set the auto-correlation length to be 0.0
  // Calculate eigenmode
  //time_t start,end;
  // Method #1: Integrate with eigenfunctions
  
   //std::time(&start);
  if(CAL_MODE){
  clock_t start = clock();
  for(int m=0;m < mode;m++){
    outputs[4][m]=0.0;
    outputs[5][m]=0.0;
    for(int i=0;i < Lx_agg; i++){
      outputs[4][m] += phi1m(m,i*Agg_Length)*height_agg[i];
      outputs[5][m] += phi2m(m,i*Agg_Length)*height_agg[i];
    }
    outputs[4][m] *= Agg_Length; // dx =  200.0
    outputs[5][m] *= Agg_Length;
   // f_phi_definition << setw(15) << outputs[4][m] << setw(15) << outputs[5][m] << "\n";
  }
  // std::time(&end);
  clock_t end = clock();
  
  // cout << "Calculating mode by definition takes " << difftime(end,start) << "\n";
  cout << "Calculating mode by definition takes " << end-start << " clocks\n";
  }
 // // Method #2: Using FFT
 // // std::time(&start);

 // // Calculate mode
 // if(CAL_MODE){
 //   start = clock();
 //   gsl_fft_real_transform(height_agg,1,Lx_agg,real,work);

 //   // The length of outputs[4][] and outputs[5][] are mode
 //   outputs[4][0] = 0.0;
 //   outputs[5][0] = SQRT_ONE_OVER_LX_AGG*height_agg[0];
 //   int m;
 //   for(m=1;m<mode-1;m++){
 //     outputs[4][m] = -SQRT_TWO_OVER_LX_AGG*height_agg[2*m];
 //     outputs[5][m] = SQRT_TWO_OVER_LX_AGG*height_agg[2*m-1];
 //   }
 //   // m = mode-1
 //   if(m==Lx_agg/2){// Using all modes, i.e., mode==Lx_agg/2+1, thus m==mode-1 ==Lx_agg/2
 //     if(Lx_agg%2==0){ // The length of input array is even
	//outputs[4][m] = 0.0;
	//outputs[5][m] = SQRT_TWO_OVER_LX_AGG*height_agg[2*m-1];
 //     }
 //     else{ // The length of input array is odd
	//outputs[4][m] = -SQRT_TWO_OVER_LX_AGG*height_agg[2*m];
	//outputs[5][m] = SQRT_TWO_OVER_LX_AGG*height_agg[2*m-1];
 //     }
 //   }
 //   else{// Use part of the mode
 //     outputs[4][m] = -SQRT_TWO_OVER_LX_AGG*height_agg[2*m];
 //     outputs[5][m] = SQRT_TWO_OVER_LX_AGG*height_agg[2*m-1];    
 //   }
  //  //for(m=0;m<mode;m++){
  //  //  f_phi_FFT << setw(15) << outputs[4][m] << setw(15) << outputs[5][m] << "\n";
  //  //}
  //  //std::time(&end);
  //  end = clock();
  //  // cout << "Calculating mode by FFT              " << difftime(end,start) << "\n";
  //  cout << "Calculating mode by FFT              " << end-start << " clocks\n";
  //}
  
  //  f_phi_definition.close();
  //  f_phi_FFT.close();

  return;
}
// --------------------------------------------------
double plant_1DSOS_RDM_input_1Dsin::getM2(){
  double Agg_Length_Square = Agg_Length*Agg_Length;
  double m2 = ((height_agg[0]-height_agg[Lx_agg-1])*(height_agg[0]-height_agg[Lx_agg-1]))/Agg_Length_Square;
  for(int i=1;i<Lx_agg;i++){
    m2+= ((height_agg[i]-height_agg[i-1])*(height_agg[i]-height_agg[i-1])/Agg_Length_Square);
  }
  m2 = m2/Lx_agg;
  return m2;
}
// --------------------------------------------------
double plant_1DSOS_RDM_input_1Dsin::getR2(){
 double mean_height = getH();
 double r2 = 0.0;
 for(int i=0;i<Lx_agg;i++){
   double dH = height_agg[i]-mean_height;
   r2 += dH*dH;
 }
 r2 = r2/Lx_agg;
 return r2;
} 
// --------------------------------------------------
double plant_1DSOS_RDM_input_1Dsin::getH(){
  double mean_height = 0;
  for(int i=0;i<Lx_agg;i++){
    mean_height += height_agg[i];
  }
  mean_height = mean_height/Lx_agg;
  return mean_height;
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::getSnapshot(const string filename){
  ofstream log_file(filename.c_str());
  for(int x=0;x<Lx;x++){
    log_file << height[x] << " " << migration_events->GetElement(x);
    if(x!=Lx-1)
      log_file << endl;
  }
  log_file.close();
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::getAggSnapshot(const string filename){
  ofstream log_file(filename.c_str());
  for(int x_agg=0;x_agg<Lx_agg;x_agg++){
    log_file << height_agg[x_agg];
    if(x_agg!=Lx_agg-1)
      log_file << endl;
  }
  log_file.close();
}

// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::set_temperature(double T_set){
  T = T_set;
  init_migr_rate();
  init_desorp_rate();

  k_m2d = exp(-(E_D-E_S)*EV2KB/(T+0.0));
  for(int x = 0;x<Lx;x++){
    int num_higher_neigh = Num_Higher_Neighbors(x);
    migration_events->SetElement(x,MIGR_RT[num_higher_neigh]);
  }
  ads_mig_des->SetElement(1,migration_events->GetTotal());
  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::set_ads_rate(double W_set,
					       double Freq_set,
					       double A_set){
  W = W_set;
  Freq = Freq_set;
  A = A_set;

  for(int x=0;x<Lx;x++){
    double W_local = W+A*sin(2*Freq*pi*(double(x)/Lx));
    adsorption_events->SetElement(x,W_local);
  }
  ads_mig_des->SetElement(0,adsorption_events->GetTotal());
  return;
}
// --------------------------------------------------
int plant_1DSOS_RDM_input_1Dsin::Num_Higher_Neighbors(int x){
  int num_higher_neighbors = 2;
  if(height[right_neigh(x)]<height[x])
    num_higher_neighbors--;
  if(height[left_neigh(x)]<height[x])
    num_higher_neighbors--;

  return num_higher_neighbors;
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::add_remove_atom(int x,int num){
  // When we add/remove an atom, the following things need to be done
  // 1. Update height
  // 2. Update migration probability of the site
  switch(num){
    case 1:
      height[x]++;break;
    case -1:
      height[x]--;break;
  }
  migration_events->SetElement(x,MIGR_RT[Num_Higher_Neighbors(x)]);
  migration_events->SetElement(left_neigh(x),MIGR_RT[Num_Higher_Neighbors(left_neigh(x))]);
  migration_events->SetElement(right_neigh(x),MIGR_RT[Num_Higher_Neighbors(right_neigh(x))]);
  ads_mig_des->SetElement(1,migration_events->GetTotal());
  return;
}
// --------------------------------------------------


// --------------------------------------------------
// Adsorption
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::adsorption(){
  add_remove_atom(adsorption_events->PickAnEvent(),1);
  return;
}
// --------------------------------------------------


// --------------------------------------------------
// Migration
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::init_migr_rate(){
  double k_m = 2*T*KB2H;
  for(int i=0;i<NUMBER_NEIGHBOR-1;i++){
    MIGR_RT[i] = k_m*exp(-(E_S+i*E_N)*EV2KB/(T+0.0));
  }
  MIGR_RT[NUMBER_NEIGHBOR-1] = 0.0;
  return;  
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::migration(){
  // Pick a site for migration
  if(migration_events->GetTotal()==0.0){
    // cerr << "Total migration rate is zero, do nothing and return\n";
    return;
  }

  int x = migration_events->PickAnEvent();
  double n;
  int destination;
  try{
    switch(Num_Higher_Neighbors(x)){
      case 0:
	add_remove_atom(x,-1);
	// Choose left or right neighbors
	n = RandN1();
	if(n>0.5){// Migrate to left neighbor  
	  add_remove_atom(left_neigh(x),1);
	}
	else{// Migrate to right neighbor
	  add_remove_atom(right_neigh(x),1);
	}
	break;
      case 1:
	
	// Find the lower neigh and migrate to it
	destination = height[left_neigh(x)]<height[x]?left_neigh(x):right_neigh(x);
	/*if((height[left_neigh(x)]<height[x]) && (height[right_neigh(x)]>=height[x])){
	  destination = left_neigh(x);
	}
	else if((height[right_neigh(x)]<height[x]) && (height[left_neigh(x)]>=height[x])){
	  destination = right_neigh(x);
	}
	else
	  throw runtime_error("Higher_Nigh is wrong");*/
	add_remove_atom(x,-1);
	add_remove_atom(destination,1);
	break;
      default:
	// throw runtime_error("Zero probability event cannot be selected");
	cerr << "One zero probability event\n";
	return;
    }
  }
  catch(runtime_error err){
    cerr << err.what();
    system("pause");
    exit(-1);
  }

  return;
}
// --------------------------------------------------


// --------------------------------------------------
// Desorption
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::init_desorp_rate(){
  double k_m = 2*T*KB2H;
  for(int i=0;i<NUMBER_NEIGHBOR;i++){
    DESORP_RT[i] = 0.0;
  }
  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM_input_1Dsin::desorption(){
  add_remove_atom(desorption_events->PickAnEvent(),-1);
  return;
}
// --------------------------------------------------


// --------------------------------------------------
// Other functions
// --------------------------------------------------
double plant_1DSOS_RDM_input_1Dsin::OneMCTime(double r_d){
  double rd;

  rd=RandRand();
  while ((rd <= 0) || (rd >= 1))
    rd=RandRand();

  return -log(rd)/r_d;
}
// --------------------------------------------------
